##############################################################################################
### This file contains analysis code for Börger, Glenk, Meyerhoff, Rehdanz "Mitigating 
### cost vector effects in stated choice experiments using cheap talk and opt-out 
### reminders" published in JAERE
###
### The associated data file is "cost-vector-data.csv"
###
### If you have comments or questions, contact Tobias Börger (tobias.boerger@hwr-berlin.de) 
###
### May 2024
##############################################################################################  



# Load packages ----------------------------------------------------------------

packages <- c("parallel", "randtoolbox", "MASS", "matrixStats", "Rfast",
              "maxLik", "numDeriv", "dplyr", "mgcv","readxl", "stringr", 
              "actuar","matlib","foreign","Matrix","readxl", "apollo", "bgw",
              "mded","stats")

lapply(packages, require, character.only = TRUE)




# A. Sample characteristics (Table 3) ------------------------------------------

database <- read.csv("cost-vector-data.csv")
data <- database[database$task==1,]

dstat <- matrix(0, ncol=8, nrow=17)

for(s in 1:8){
  dstat[1,s]  <- 100*sum(data$female[data$split==s])/length(data$female[data$split==s])
  dstat[2,s]  <- mean(data$age[data$split==s])
  dstat[3,s]  <- 100*(table(factor(data$v_375[data$split==s], levels=c(1,2,3,4,6))) / length(data$v_375[data$split==s]))[1]
  dstat[4,s]  <- 100*(table(factor(data$v_375[data$split==s], levels=c(1,2,3,4,6))) / length(data$v_375[data$split==s]))[2]
  dstat[5,s]  <- 100*(table(factor(data$v_375[data$split==s], levels=c(1,2,3,4,6))) / length(data$v_375[data$split==s]))[3]
  dstat[6,s]  <- 100*(table(factor(data$v_375[data$split==s], levels=c(1,2,3,4,6))) / length(data$v_375[data$split==s]))[4]
  dstat[7,s]  <- 100*(table(factor(data$v_375[data$split==s], levels=c(1,2,3,4,6))) / length(data$v_375[data$split==s]))[5]
  dstat[8,s]  <- 100*(table(factor(data$v_380[data$split==s], levels=0:9)) / length(data$v_380[data$split==s]))[1]
  dstat[9,s]  <- 100*(table(factor(data$v_380[data$split==s], levels=0:9)) / length(data$v_380[data$split==s]))[2]
  dstat[10,s] <- 100*(table(factor(data$v_380[data$split==s], levels=0:9)) / length(data$v_380[data$split==s]))[3]
  dstat[11,s] <- 100*(table(factor(data$v_380[data$split==s], levels=0:9)) / length(data$v_380[data$split==s]))[4]
  dstat[12,s] <- 100*(table(factor(data$v_380[data$split==s], levels=0:9)) / length(data$v_380[data$split==s]))[5]
  dstat[13,s] <- 100*(table(factor(data$v_380[data$split==s], levels=0:9)) / length(data$v_380[data$split==s]))[6]
  dstat[14,s] <- 100*(table(factor(data$v_380[data$split==s], levels=0:9)) / length(data$v_380[data$split==s]))[7]
  dstat[15,s] <- 100*(table(factor(data$v_380[data$split==s], levels=0:9)) / length(data$v_380[data$split==s]))[8]
  dstat[16,s] <- 100*(table(factor(data$v_380[data$split==s], levels=0:9)) / length(data$v_380[data$split==s]))[9]
  dstat[17,s] <- 100*(table(factor(data$v_380[data$split==s], levels=0:9)) / length(data$v_380[data$split==s]))[10]
}
rownames(dstat) <- c("female (percent)","age (mean)","Haupt/Real","Ausbildung","Abitur","Universität","Andere",
                     "Inc: 0","Inc: <€900","Inc: €900-1299","Inc: €1300-1499","Inc: €1500-1999","Inc: €2000-2599","Inc: €2600-3599","Inc: €3600-4999","Inc: >€5000","Inc: n/a")
colnames(dstat) <- c("T1","T2","T3","T4","T5","T6","T7","T8")

### Table 3 ###
round(dstat, digits=2)

# ####



# B. share of SQ per treatment (Table 4) ---------------------------------------

database <- read.csv("cost-vector-data.csv")

choices <- matrix(database$choi, nrow=nrow(database)/8, ncol=8, byrow=TRUE)
nsq <- rowCounts(choices, value=1)
choices <- as.data.frame(cbind(newid=database$newid[seq(1, length(database$newid), 8)], split=database$split[seq(1, length(database$split), 8)], choices, nsq=nsq))

table1 <- matrix(c(table(factor(choices$nsq[choices$split==1], levels=0:8)),
                   table(factor(choices$nsq[choices$split==2], levels=0:8)),
                   table(factor(choices$nsq[choices$split==3], levels=0:8)),
                   table(factor(choices$nsq[choices$split==4], levels=0:8)),
                   table(factor(choices$nsq[choices$split==5], levels=0:8)),
                   table(factor(choices$nsq[choices$split==6], levels=0:8)),
                   table(factor(choices$nsq[choices$split==7], levels=0:8)),
                   table(factor(choices$nsq[choices$split==8], levels=0:8))), nrow=9, ncol=8, byrow=FALSE)  
table1 <- rbind(table1, colSums(table1))
colnames(table1) <- rep(c("Cost300","Cost500","Cost1000","Cost1800"),2)
rownames(table1) <- c("Never","1","2","3","4","5","6","7","Always","Sum")
table1

table2 <- table1
for(i in 1:ncol(table2)){
  table2[,i] <- table1[,i] / table1[10,i]
}

### Table 4 ###
round(table2*100, digits=2)
# ####



# C. cumulative bid acceptance curves (Table A.8, A.9, Figure 3, 5, A.1) -------

database <- read.csv("cost-vector-data.csv")

cum.accept.bid <- matrix(0, nrow=8, ncol=9)
for(s in 1:8){
  cost.levs <- c(0,sort(unique(database$cost1[database$split==s])))
  for(l in 9:1){
    choices_temp1 <- ifelse(database$choice1[database$split==s]==1 & database$cost1[database$split==s]>=cost.levs[l],1,0)
    choices_temp2 <- ifelse(database$choice2[database$split==s]==1 & database$cost2[database$split==s]>=cost.levs[l],1,0)
    cum.accept.bid[s,l] <- (sum(choices_temp1)+sum(choices_temp2)) / nrow(database[database$split==s,]) 
  }
}
cum.accept.bid[,1] <- 1
colnames(cum.accept.bid) <- c("zero","lev1","lev2","lev3","lev4","lev5","lev6","lev7","lev8")
rownames(cum.accept.bid) <- c("t1","t2","t3","t4","t5","t6","t7","t8")

### Table A.8 ###
round(cum.accept.bid*100, digits=2)


### Figure 3 ###
col <- c("maroon","darkorchid3","blue3","forestgreen")

par(mfrow=c(2,1), mar=c(3.15,4,0.5,0.2))
plot(0,0, col="white", xlim=c(0,1800), ylim=c(0,100), ylab="Cummulative bid acceptance (in percent)", las=1)
abline(h=seq(0,100,10), col="grey")
points(c(0,sort(unique(database$cost1[database$split==4]))),cum.accept.bid[4,]*100, lwd=3, type="b", pch=16, col=col[4])
points(c(0,sort(unique(database$cost1[database$split==3]))),cum.accept.bid[3,]*100, lwd=3, type="b", pch=15, col=col[3])
points(c(0,sort(unique(database$cost1[database$split==2]))),cum.accept.bid[2,]*100, lwd=3, type="b", pch=1, col=col[2])
points(c(0,sort(unique(database$cost1[database$split==1]))),cum.accept.bid[1,]*100, lwd=3, type="b", pch=17, col=col[1])
title(xlab="Cost vector element (in EUR)", line=2.1, cex.lab=1)
legend(x="topright", legend=c("Cost300","Cost500","Cost1000","Cost1800"), cex=1.5, col=col, lwd=rep(3,4), pch=c(17,1,15,16), border="white")

plot(0,0, col="white", xlim=c(0,1800), ylim=c(0,100), ylab="Cummulative bid acceptance (in percent)", las=1)
abline(h=seq(0,100,10), col="grey")
points(c(0,sort(unique(database$cost1[database$split==8]))),cum.accept.bid[8,]*100, lwd=3, type="b", pch=16, col=col[4])
points(c(0,sort(unique(database$cost1[database$split==7]))),cum.accept.bid[7,]*100, lwd=3, type="b", pch=15, col=col[3])
points(c(0,sort(unique(database$cost1[database$split==6]))),cum.accept.bid[6,]*100, lwd=3, type="b", pch=1, col=col[2])
points(c(0,sort(unique(database$cost1[database$split==5]))),cum.accept.bid[5,]*100, lwd=3, type="b", pch=17, col=col[1])
title(xlab="Cost vector element (in EUR)", line=2.1, cex.lab=1)
legend(x="topright", legend=c("Cost300","Cost500","Cost1000","Cost1800"), cex=1.5, col=col, lwd=rep(3,4), pch=c(17,1,15,16), border="white")



### Figure A.1 ###
col <- c("maroon","darkorchid3","blue3","forestgreen")

par(mfrow=c(2,2), mar=c(4.5,4.3,0.5,0.2))
plot(0, col="white", las=1, xlim=c(0,300), ylim=c(0,100), xlab="Cost amount (in EUR)", ylab="Cum. bid acceptance (in percent)", cex.lab=1.4, cex.axis=1.4)
abline(h=seq(0,100,10), col="grey")
points(c(0,sort(unique(database$cost1[database$split==1]))),cum.accept.bid[1,]*100, lwd=3, type="b", pch=16, col=col[1])
points(c(0,sort(unique(database$cost1[database$split==5]))),cum.accept.bid[5,]*100, lwd=3, type="b", pch=4,  col=col[1])
legend(x="topright", legend=c("without device","with CT&OOR device"), cex=1.4, col=col[c(1,1)], lwd=rep(2,4), pch=c(16,4), border="white")

plot(0, col="white", las=1, xlim=c(0,500), ylim=c(0,100), xlab="Cost amount (in EUR)", ylab="Cum. bid acceptance (in percent)", cex.lab=1.4, cex.axis=1.4)
abline(h=seq(0,100,10), col="grey")
points(c(0,sort(unique(database$cost1[database$split==2]))),cum.accept.bid[2,]*100, lwd=3, type="b", pch=16, col=col[2])
points(c(0,sort(unique(database$cost1[database$split==6]))),cum.accept.bid[6,]*100, lwd=3, type="b", pch=4, col=col[2])
legend(x="topright", legend=c("without device","with CT&OOR device"), cex=1.4, col=col[c(2,2)], lwd=rep(2,4), pch=c(16,4), border="white")

plot(0, col="white", las=1, xlim=c(0,1000), ylim=c(0,100), xlab="Cost amount (in EUR)", ylab="Cum. bid acceptance (in percent)", cex.lab=1.4, cex.axis=1.4)
abline(h=seq(0,100,10), col="grey")
points(c(0,sort(unique(database$cost1[database$split==3]))),cum.accept.bid[3,]*100, lwd=3, type="b", pch=16, col=col[3])
points(c(0,sort(unique(database$cost1[database$split==7]))),cum.accept.bid[7,]*100, lwd=3, type="b", pch=4,  col=col[3])
legend(x="topright", legend=c("without device","with CT&OOR device"), cex=1.4, col=col[c(3,3)], lwd=rep(2,4), pch=c(16,4), border="white")

plot(0, col="white", las=1, xlim=c(0,1800), ylim=c(0,100), xlab="Cost amount (in EUR)", ylab="Cum. bid acceptance (in percent)", cex.lab=1.4, cex.axis=1.4)
abline(h=seq(0,100,10), col="grey")
points(c(0,sort(unique(database$cost1[database$split==4]))),cum.accept.bid[4,]*100, lwd=3, type="b", pch=16, col=col[4])
points(c(0,sort(unique(database$cost1[database$split==8]))),cum.accept.bid[8,]*100, lwd=3, type="b", pch=4, col=col[4])
legend(x="topright", legend=c("without device","with CT&OOR device"), cex=1.4, col=col[c(4,4)], lwd=rep(2,4), pch=c(16,4), border="white")




# Bid acceptance for high- and low-income respondents

database <- read.csv("cost-vector-data.csv")

inc.gr <- list(c(0,1,2,3,4) ,c(5,6,7,8))

cum.accept.bid_temp <- matrix(1:(4*8), nrow=4, ncol=8)
accept.bid.count_temp <- all.choices.count_temp <- matrix(7, nrow=4, ncol=8)

for(k in 1:2){
  # source(file.path("data","setup-data.R"))
  database <- read.csv("cost-vector-data.csv")
  database <- database[database$v_380 %in% inc.gr[[k]], ]
  for(s in 1:2){
    for(l in 8:1){
      choices_temp1 <- ifelse(database$choice1[database$CheapT==s-1]==1 & database$costlev1[database$CheapT==s-1]>=l,1,0)
      choices_temp2 <- ifelse(database$choice2[database$CheapT==s-1]==1 & database$costlev2[database$CheapT==s-1]>=l,1,0)
      cum.accept.bid_temp[(k-1)*2+s,l] <- (sum(choices_temp1)+sum(choices_temp2)) / nrow(database[database$CheapT==s-1,]) 
      # Store total # of non-sq choice and total # of choices for hypothesis testing (Chi2)
      accept.bid.count_temp[(k-1)*2+s,l] <- sum(choices_temp1) + sum(choices_temp2)  
      all.choices.count_temp[(k-1)*2+s,l] <- nrow(database[database$CheapT==s-1,])
    }
  }
}
cum.accept.bid <- cum.accept.bid_temp[c(1,3,2,4),]
cum.accept.bid <- cbind(rep(1, nrow(cum.accept.bid)), cum.accept.bid)
colnames(cum.accept.bid) <- c("zero","lev1","lev2","lev3","lev4","lev5","lev6","lev7","lev8")
rownames(cum.accept.bid) <- c("no device-low.inc","no device-high.inc","CT&OOR-low.inc","CT&OOR-high.inc")
round(cum.accept.bid, 3)

chi2.test.pvals <- matrix(3, nrow=2, ncol=8)  # Chi-squared tests

accept.bid.count  <- accept.bid.count_temp[c(1,3,2,4),]
all.choices.count <- all.choices.count_temp[c(1,3,2,4),]
colnames(accept.bid.count) <- colnames(all.choices.count) <- colnames(chi2.test.pvals) <- c("lev1","lev2","lev3","lev4","lev5","lev6","lev7","lev8")
rownames(accept.bid.count) <- rownames(all.choices.count) <- c("no device-low.inc","no device-high.inc","CT&OOR-low.inc","CT&OOR-high.inc")
rownames(chi2.test.pvals) <- c("Comparison (no device)", "Comparison CT&OOR)")

for(r in 1:nrow(chi2.test.pvals)){
  for(c in 1:ncol(chi2.test.pvals)){
    chi2.test.pvals[r,c] <- chisq.test(rbind(c(accept.bid.count[2*r-1,c],all.choices.count[2*r-1,c]-accept.bid.count[2*r-1,c]),
                                             c(accept.bid.count[2*r,  c],all.choices.count[2*r,  c]-accept.bid.count[2*r,  c])))$p.value
  }
}

### Table A.9 ###
round(cum.accept.bid, 3)
round(chi2.test.pvals, 3)



### Figure 5 ###
col <- rep("black", 2) 
lwd <- 2

par(mfrow=c(2,1), mar=c(2.5,4,0.5,0.2))
plot(cum.accept.bid[1,-1]*100, type="b", col=col[1], pch=16, ylim=c(0,90), lwd=lwd, las=1, xlab="Bid level"
     , ylab="Cummulative bid acceptance (in %)")
lines(cum.accept.bid[2,-1]*100, type="b", col=col[2], pch=16, lwd=lwd, lty=5)
text(x=5, y=85, "A. no CT&OOR device (T1 to T4)", cex=1.6)
legend(x="bottomleft", legend=c("Low income","High income"), cex=1.3, col=col, lwd=rep(lwd,2), lty=c(1,5), border="white")

plot(cum.accept.bid[3,-1]*100, type="b", col=col[1], pch=16, ylim=c(0,90), lwd=lwd, las=1, xlab="Bid level"
     , ylab="Cummulative bid acceptance (in %)")
lines(cum.accept.bid[4,-1]*100, type="b", col=col[2], pch=16, lwd=lwd, lty=5)
text(x=5, y=85, "B. CT&OOR device included (T5 to T8)", cex=1.6)
legend(x="bottomleft", legend=c("Low income","High income"), cex=1.3, col=col, lwd=rep(lwd,2), lty=c(1,5), border="white")
# ####



# D. MXL models in WTP space (Tables A.1, A.2) ---------------------------------

rm(choice0,choice1,choice2)

mods <- list()

for(s in 1:8){
  
  #### 1. LOAD LIBRARY AND DEFINE CORE SETTINGS  ####
  ### Initialise code
  apollo_initialise()
  
  ### Set core controls
  apollo_control = list(
    modelName  ="MXL-wtpspace",
    modelDescr ="MXL model",
    indivID    ="newid",
    mixing     = TRUE,
    nCores     = 10,
    panelData  = TRUE
  )
  # ####
  
  #### 2. LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
  database <- read.csv("cost-vector-data.csv")
  database <- database[database$split==s, ]
  database$cost1 <- database$cost1/100
  database$cost2 <- database$cost2/100
  # ####
  
  #### 3. DEFINE MODEL PARAMETERS   ####
  if(s==1){
    apollo_beta <- c(-2.65263208635266, 0.405952982477775, 0.767774545705003, 1.01921575094023, 0.481923289717677, 0.602285530453235, 1.65506194441807, 1.4390581237503, 0.0559832335228737, 0.792651448852516, 1.8392982954345, 1.54978680287548, 0.541880417919293, 4.6259552827646, -0.204982843165817, 0.830837427330742, 0.618580655653492, -0.53364627573815, -0.303995207116986, -0.590946738051033, -0.570532101434146, -0.34064745867358, 0.789734899133393, 0.950584120935758, -0.731980573350487, -1.47522689820307)
  } else if(s==2){
    apollo_beta <- c(-2.77620745910405, 0.437888852605115, 0.463020857664999, 1.22770051207973, 0.668659866094868, 0.556413734747439, 1.91169721216585, 1.96679259026373, 0.0591698402016176, 1.96777012775607, 2.9791652851142, 2.14928650342964, 0.12069851483603, 4.11311386005255, 0.503029170855143, 0.760677742345328, -0.992672791126551, -0.540764993850944, 0.871586991785066, -0.514162740733948, 0.985324925171082, -0.457047181819632, -0.320201805318167, -1.77418849503464, -1.28358453655197, 1.74244620390512)
  } else if(s==3){
    apollo_beta <- c(-3.17086552781122, -1.09923218309314, -0.883053341273167, 2.28392629590673, 1.81443584496903, 0.653122253979749, 3.71358155186659, 3.44317138899486, -0.128759772525961, 2.27474336622064, 4.58235672206017, 3.94761029863487, -0.86160491159508, -4.14894826832152, 1.0931804886051, -1.35703563698518, 2.32676297154353, -1.83746169078137, -1.22676042148575, -0.667284759935864, -2.19058349781554, -1.24612342779362, 3.42776409033132, 3.24525872904867, 3.03043393532906, 0.463282725151401)
  } else if(s==4){
    apollo_beta <- c(-3.29239539149645, -4.84778498375857, -3.33492449936475, 2.37140387498517, 3.92232687501133, 1.38290399709659, 11.2649913921298, 6.12084537993523, -0.828452454428715, 3.42538000947281, 6.72748072978897, 5.72971608287544, -1.75963736874492, -4.16347342843604, 3.74001430137861, -6.25848800704848, 3.63233535743214, -4.0235253970223, 8.31210292885724, 6.8014903116618, -5.46266228193461, -0.121658495448265, 8.20650157725489, -2.95408549361164, 7.22317929844729, 0.702537151596151)
  } else if(s==5){
    apollo_beta <- c(-1.76263354824628, 0.13188107823322, 0.216936783160683, 0.360594414146951, 0.107918366571101, -0.0419166764115528, 0.239767777656736, 0.35371625389067, -0.0258471038697665, 0.0559523843834879, 0.222235912687907, 0.26012304301049, 0.431819504473352, 3.28187134867306, 0.199437915781725, 0.0916200897328642, -0.237501254412381, 0.0433054643812685, 0.0131649552979268, -0.13149469991346, 0.0817712764797097, 0.02156620928175, 0.0580608987647965, 0.0711623628687085, 0.0532933596340941, 2.43709954391831)
  } else if(s==6){
    apollo_beta <- c(-1.46099323203738, 0.186091749433099, 0.131393336961229, 0.426685687729675, 0.0579722907943861, -0.0825950572377967, 0.306523312444835, 0.460921025804731, -0.0368195238472204, 0.0956691219101148, 0.304511263978093, 0.347202976833434, 0.330227868000204, 3.31835804454631, 0.127627973132148, 0.00975738615672001, 0.294598714453172, -0.109516982760388, -0.0313224751408578, 0.0161118097663892, -0.143495528423689, 0.000485020011957919, 0.0560044370240051, -0.101309176807532, 0.156784285079606, 2.20645150480676)
  } else if(s==7){
    apollo_beta <- c(-1.50903023769313, 0.138753071104559, 0.196038554488902, 0.369412850969157, 0.0430352401119843, -0.0145149651947564, 0.257397652907333, 0.360175336044234, -0.0711141171029675, 0.137176399299319, 0.204874908309153, 0.283086605589281, 0.102615714295952, 3.52470265999877, 0.00230116124654033, 0.00551379316520591, 0.119709944910345, 0.0356261828343448, 0.0122320859861843, -0.163933695234377, 0.140664983504661, -0.0225649115511815, 3.75161508438519e-05, 0.0886229356346607, 0.181688946704, 2.36792594780692)
  } else if(s==8){
    apollo_beta <- c(-1.3926896346499, 0.473935862333245, 0.295392693978071, 0.342232669471805, 0.121735493439782, 0.000165828601352884, 0.420352020303707, 0.501096020298621, -0.0127151230937205, -0.0786039443055845, 0.200061324676577, 0.19572129747957, -0.197894172514738, 2.99249352565368, 0.196502693720949, 0.23451824828637, -0.0131148157628762, 0.00311281758348218, -0.0105589871137515, 0.171762379243562, 0.0102937508978971, -0.0106853548809317, 0.083708916308781, -0.0143450088741923, 0.133683886871985, 2.69577388130183)
  }
  
  names(apollo_beta) <- c("mu_none","mu_clar1","mu_clar2","mu_clar3","mu_fish","mu_bio1","mu_bio2","mu_bio3",
                          "mu_coast","mu_lit1","mu_lit2","mu_lit3","mu_cost",
                          "sigma_none","sigma_clar1","sigma_clar2","sigma_clar3","sigma_fish","sigma_bio1","sigma_bio2","sigma_bio3",
                          "sigma_coast","sigma_lit1","sigma_lit2","sigma_lit3","sigma_cost")
  apollo_fixed = c()  
  # ####
  
  #### 4. DEFINE RANDOM COMPONENTS                                    ####
  ### Set parameters for generating draws
  apollo_draws = list(
    interDrawsType = "sobol",
    interNDraws    = 1000,
    interNormDraws = c("draws_none", "draws_clar_1", "draws_clar_2", "draws_clar_3", "draws_fish",
                       "draws_bio_1", "draws_bio_2", "draws_bio_3", "draws_coast",
                       "draws_lit_1", "draws_lit_2", "draws_lit_3", "draws_cost")
  )
  
  ### Create random parameters
  apollo_randCoeff = function(apollo_beta, apollo_inputs){
    randcoeff = list()
    
    randcoeff[["b_cost"]]  = -exp(mu_cost + sigma_cost*draws_cost)
    randcoeff[["b_none"]]  = mu_none  + sigma_none*draws_none    
    randcoeff[["b_clar1"]] = mu_clar1 + sigma_clar1*draws_clar_1 
    randcoeff[["b_clar2"]] = mu_clar2 + sigma_clar2*draws_clar_2 
    randcoeff[["b_clar3"]] = mu_clar3 + sigma_clar3*draws_clar_3 
    randcoeff[["b_fish"]]  = mu_fish  + sigma_fish*draws_fish    
    randcoeff[["b_bio1"]]  = mu_bio1  + sigma_bio1*draws_bio_1   
    randcoeff[["b_bio2"]]  = mu_bio2  + sigma_bio2*draws_bio_2   
    randcoeff[["b_bio3"]]  = mu_bio3  + sigma_bio3*draws_bio_3   
    randcoeff[["b_coast"]] = mu_coast + sigma_coast*draws_coast  
    randcoeff[["b_lit1"]]  = mu_lit1  + sigma_lit1*draws_lit_1   
    randcoeff[["b_lit2"]]  = mu_lit2  + sigma_lit2*draws_lit_2   
    randcoeff[["b_lit3"]]  = mu_lit3  + sigma_lit3*draws_lit_3   
    
    return(randcoeff)
  }
  # ####
  
  #### 5. GROUP AND VALIDATE INPUTS                                   ####
  apollo_inputs = apollo_validateInputs()
  # ####
  
  #### 6. DEFINE MODEL AND LIKELIHOOD FUNCTION                        ####
  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
    
    ### Attach inputs and detach after function exit
    apollo_attach(apollo_beta, apollo_inputs)
    on.exit(apollo_detach(apollo_beta, apollo_inputs))
    
    ### Create list of probabilities P
    P = list()
    ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
    V = list()
    V[['sq']]   =  b_none + b_cost*(cost0 -( b_clar1*clar1_0 + b_clar2*clar2_0 + b_clar3*clar3_0 + b_fish*fish0 + b_bio1*bio1_0 + b_bio2*bio2_0 + b_bio3*bio3_0 + b_coast*coast0 + b_lit1*lit1_0 + b_lit2*lit2_0 + b_lit3*lit3_0))
    V[['alt2']] =           b_cost*(cost1 -( b_clar1*clar1_1 + b_clar2*clar2_1 + b_clar3*clar3_1 + b_fish*fish1 + b_bio1*bio1_1 + b_bio2*bio2_1 + b_bio3*bio3_1 + b_coast*coast1 + b_lit1*lit1_1 + b_lit2*lit2_1 + b_lit3*lit3_1))
    V[['alt3']] =           b_cost*(cost2 -( b_clar1*clar1_2 + b_clar2*clar2_2 + b_clar3*clar3_2 + b_fish*fish2 + b_bio1*bio1_2 + b_bio2*bio2_2 + b_bio3*bio3_2 + b_coast*coast2 + b_lit1*lit1_2 + b_lit2*lit2_2 + b_lit3*lit3_2))
    
    ### Define settings for MNL model component
    mnl_settings = list(
      alternatives = c(sq=1, alt2=2, alt3=3),
      avail        = 1,
      choiceVar    = choi,
      V            = V
    )
    
    ### Compute probabilities using MNL model
    P[['model']] = apollo_mnl(mnl_settings, functionality)
    
    ### Take product across observation for same individual     
    P = apollo_panelProd(P, apollo_inputs, functionality)
    
    ### Average across inter-individual draws
    P = apollo_avgInterDraws(P, apollo_inputs, functionality)
    
    ### Prepare and return outputs of function
    P = apollo_prepareProb(P, apollo_inputs, functionality)
    return(P)
  }
  # ####
  
  #### 7. MODEL ESTIMATION                                            ####
  estimate_settings = list(estimationRoutine="bgw", hessianRoutine="numDeriv")
  mods[[s]] = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, estimate_settings)
  # #####
}

### (values for) Table A.1 ###
summary(mods[[1]])
summary(mods[[2]])
summary(mods[[3]])
summary(mods[[4]])

### (values for) Table A.2 ###
summary(mods[[5]])
summary(mods[[6]])
summary(mods[[7]])
summary(mods[[8]])



# ALTERNATIVELY: Read in saved model estimates and continue postestimation in Secion E

mods <- readRDS("MXL-AllTreat.rds")

### (values for) Table A.1 ###
summary(mods[[1]])
summary(mods[[2]])
summary(mods[[3]])
summary(mods[[4]])

### (values for) Table A.2 ###
summary(mods[[5]])
summary(mods[[6]])
summary(mods[[7]])
summary(mods[[8]])

# ####



# E. WTP (Table 5, A.3, A.4, and Figure 4) -------------------------------------

# compute WTP and simulate confidence intervals 
ndraws <- 10000
wtp.md <- wtp.lb <- wtp.ub <- wtp.sd <- matrix(0, nrow=8, ncol=12)
draws.mods <- list()

for(s in 1:8){
  model_temp <- mods[[s]]
  wtp.md[s,] <- model_temp$estimate[1:12] * 100
  draws <- mvrnorm(ndraws, mu=model_temp$estimate, Sigma=model_temp$robvarcov)
  draws.mods[[s]] <- draws
  wtp_sim <- matrix(0, nrow=ndraws, ncol=ncol(wtp.md))
  for(i in 1:ndraws){
    wtp_sim[i,] <- draws[i,1:12] * 100
  }
  for(a in 1:ncol(wtp.md)){
    wtp.lb[s,a] <- quantile(wtp_sim[,a], probs=0.025)
    wtp.ub[s,a] <- quantile(wtp_sim[,a], probs=0.975)
    wtp.sd[s,a] <- sd(wtp_sim[,a])*100
  }
}
names.wtp <- c("None","clar1","clar2","clar3","fish","bio1","bio2","bio3","coast","lit1","lit2","lit3")
names.treat <- c("T1","T2","T3","T4","T5","T6","T7","T8")
colnames(wtp.md) <- colnames(wtp.lb) <- colnames(wtp.ub) <- names.wtp
rownames(wtp.md) <- rownames(wtp.lb) <- rownames(wtp.ub) <- names.treat

### (figures for) Tables 5, A.3 and A.4 ###
round(wtp.md, digits=2)        
round(wtp.lb, digits=2)        
round(wtp.ub, digits=2) 

round(wtp.sd, digits=2)
round(abs(wtp.sd/wtp.md), digits=2)


# plot relative changes in mWTP (Figure 4)
toplot1 <- toplot2 <- ci.low1 <- ci.high1 <- ci.low2 <- ci.high2 <- matrix(0, nrow=11, ncol=4)
for(i in 1:11){
  for(s in 1:4){
    toplot1[i,s] <- wtp.md[s,i+1] / wtp.md[1,i+1]
    ci.low1[i,s] <- quantile(draws.mods[[s]][,i+1] / draws.mods[[1]][,i+1], probs=0.025)
    ci.high1[i,s] <- quantile(draws.mods[[s]][,i+1] / draws.mods[[1]][,i+1], probs=0.975)
    toplot2[i,s] <- wtp.md[s+4,i+1] / wtp.md[5,i+1]
    ci.low2[i,s] <- quantile(draws.mods[[s+4]][,i+1] / draws.mods[[5]][,i+1], probs=0.025)
    ci.high2[i,s] <- quantile(draws.mods[[s+4]][,i+1] / draws.mods[[5]][,i+1], probs=0.975)    
  }
}

toplot1[1:2,3:4] <- toplot1[5,3:4] <- toplot1[8,3:4] <- NA
toplot2[c(5,8),] <- toplot2[9,4] <- NA
ci.low1[1:2,3:4] <- ci.low1[5,3:4] <- ci.low1[8,3:4] <- NA
ci.low2[c(5,8),] <- ci.low2[9,4] <- NA
ci.high1[1:2,3:4] <- ci.high1[5,3:4] <- ci.high1[8,3:4] <- NA
ci.high2[c(5,8),] <- ci.high2[9,4] <- NA


### Figure 4 ###
col <- rep("black",2)
lwd <- 3
cex <- 1.6

par(mfrow=c(2,2), mar=c(2.4,4.2,1.8,0.2), cex.lab=1.3)
plot( 1:4, toplot1[1,]*100, type="b", col=col[1], lwd=lwd, pch=16, cex=cex, xaxt="n", yaxt="n", xlim=c(0.8,12.2), 
      ylim=c(0,820), main="Clarity", las=2, ylab="Marg. WTP rel. to treatment 'Cost300'")
lines(1:4,  toplot2[1,]*100, type="b", col=col[2], lwd=lwd, pch=1, cex=cex)
lines(5:8,  toplot1[2,]*100, type="b", col=col[1], lwd=lwd, pch=16, cex=cex)
lines(5:8,  toplot2[2,]*100, type="b", col=col[2], lwd=lwd, pch=1, cex=cex)
lines(9:12, toplot1[3,]*100, type="b", col=col[1], lwd=lwd, pch=16, cex=cex)
lines(9:12, toplot2[3,]*100, type="b", col=col[2], lwd=lwd, pch=1, cex=cex)
axis(side=2,at=seq(0,800,100), las=1)
axis(side=1,at=c(2.5,6.5,10.5), c("rather turbid","rather clear","clear"), cex.axis=1.4, tick=F)
abline(h=100, v=c(4.5,8.5))

plot( 1:4, toplot1[5,]*100, type="b", col=col[1], lwd=lwd, pch=16, cex=cex, xaxt="n", yaxt="n", xlim=c(0.8,12.2), 
      ylim=c(0,820), main="State of biodiversity", las=2, ylab="Marg. WTP rel. to treatment 'Cost300'")
lines(1:4,  toplot2[5,]*100, type="b", col=col[2], lwd=lwd, pch=1, cex=cex)
lines(5:8,  toplot1[6,]*100, type="b", col=col[1], lwd=lwd, pch=16, cex=cex)
lines(5:8,  toplot2[6,]*100, type="b", col=col[2], lwd=lwd, pch=1, cex=cex)
lines(9:12, toplot1[7,]*100, type="b", col=col[1], lwd=lwd, pch=16, cex=cex)
lines(9:12, toplot2[7,]*100, type="b", col=col[2], lwd=lwd, pch=1, cex=cex)
axis(side=2,at=seq(0,800,100), las=1)
axis(side=1,at=c(2.5,6.5,10.5), c("rather bad","rather good","good"), cex.axis=1.4, tick=F)
abline(h=100, v=c(4.5,8.5))

plot( 1:4, toplot1[4,]*100, type="b", col=col[1], lwd=lwd, pch=16, cex=cex, xaxt="n", yaxt="n", xlim=c(0.8,8.2), 
      ylim=c(0,820), main="Fish stocks / Coastal protection", las=2, ylab="Marg. WTP rel. to treatment 'Cost300'")
lines(1:4,  toplot2[4,]*100, type="b", col=col[2], lwd=lwd, pch=1, cex=cex)
lines(5:8,  toplot1[8,]*100, type="b", col=col[1], lwd=lwd, pch=16, cex=cex)
lines(5:8,  toplot2[8,]*100, type="b", col=col[2], lwd=lwd, pch=1, cex=cex)
axis(side=2,at=seq(0,800,100), las=1)
axis(side=1,at=c(2.5,6.5), c("stable fish stocks","mainly low impact"), cex.axis=1.4, tick=F)
abline(h=100, v=4.5)

plot( 1:4, toplot1[9,]*100, type="b", col=col[1], lwd=lwd, pch=16, cex=cex, xaxt="n", yaxt="n", xlim=c(0.8,12.2), 
      ylim=c(0,820), main="Amount of litter", las=2, ylab="Marg. WTP rel. to treatment 'Cost300'")
lines(1:4,  toplot2[9,]*100,  type="b", col=col[2], lwd=lwd, pch=1, cex=cex)
lines(5:8,  toplot1[10,]*100, type="b", col=col[1], lwd=lwd, pch=16, cex=cex)
lines(5:8,  toplot2[10,]*100, type="b", col=col[2], lwd=lwd, pch=1, cex=cex)
lines(9:12, toplot1[11,]*100, type="b", col=col[1], lwd=lwd, pch=16, cex=cex)
lines(9:12, toplot2[11,]*100, type="b", col=col[2], lwd=lwd, pch=1, cex=cex)
axis(side=2,at=seq(0,800,100), las=1)
axis(side=1,at=c(2.5,6.5,10.5), c("rather bad","rather good","good"), cex.axis=1.4, tick=F)
abline(h=100, v=c(4.5,8.5))
# ####



# F. WTP post-estimation (Tables A.5, A.6) ------------------------------------- 

# tests differences in WTP across treatments using Poe et al. (2005) 

ndraws <- 10000
wtp.md <- matrix(0, nrow=8, ncol=11)
wtp_sim <- list()
for(s in 1:8){
  model_temp <- mods[[s]]
  wtp.md[s,] <- model_temp$estimate[2:12] * 100
  set.seed(123476) ; draws <- mvrnorm(ndraws, mu=model_temp$estimate, Sigma=model_temp$robvarcov)
  wtp_sim[[s]] <- draws[,2:12] * 100
}

poe.res <- list()

for(a in 1:11){
  att <- a
  poe_dif <- matrix(9, nrow=8, ncol=8)
  poe_res <- matrix(9, nrow=8, ncol=8)
  
  for(i in 1:8){
    for(j in 1:8){
      if (i==j) {
        poe_res[i,j] <- NA
      } else {
        if(mean(wtp_sim[[i]][,att]) > mean(wtp_sim[[j]][,att])) {
          dist1 <- wtp_sim[[i]][,att]
          dist2 <- wtp_sim[[j]][,att] 
          test_temp <- mded(distr1=dist1, distr2=dist2, detail=FALSE)
          poe_res[i,j] <- test_temp$stat
          poe_dif[i,j] <- test_temp$means[1] - test_temp$means[2]
        } else if (mean(wtp_sim[[j]][,att]) > mean(wtp_sim[[i]][,att])) {
          dist2 <- wtp_sim[[i]][,att]
          dist1 <- wtp_sim[[j]][,att] 
          test_temp <- mded(distr1=dist1, distr2=dist2, detail=FALSE)
          poe_res[i,j] <- test_temp$stat
          poe_dif[i,j] <- test_temp$means[1] - test_temp$means[2]
        }
        cat("\r",(i-1)*8+j,"of",64,"possible tests completed, for attribute",a,"of 11")
      }
    }
  }
  poe.res[[a]] <- poe_res
}


# Benjamini-Hochberg (1995) correction of p-values for treatment arms separately

label1 <- rep(c("clar1","clar2","clar3","fish","bio1","bio2","bio3","coast","lit1","lit2","lit3"), each=6) 
label2_1 <- rep(c("T1 vs. T2","T1 vs. T3","T1 vs. T4","T2 vs. T3","T2 vs. T4","T3 vs. T4"), 11)  
label2_2 <- rep(c("T5 vs. T6","T5 vs. T7","T5 vs. T8","T6 vs. T7","T6 vs. T8","T7 vs. T8"), 11)  

# extract p-value for treatment arm 1
pval1_mat <- pval2_mat <- matrix(9, nrow=11, ncol=6)

for(i in 1:11){
  pval1_mat[i,] <- as.vector(poe.res[[i]][1:4,1:4][lower.tri(poe.res[[i]][1:4,1:4])==TRUE])         
  pval2_mat[i,] <- as.vector(poe.res[[i]][5:8,5:8][lower.tri(poe.res[[i]][5:8,5:8])==TRUE])         
}

pval1 <- as.vector(t(pval1_mat))
pval2 <- as.vector(t(pval2_mat))

ranks1 <- rank(pval1, ties.method="last")
ranks2 <- rank(pval2, ties.method="last")

p_m_over_k1 <- pval1*length(pval1)/ranks1
p_m_over_k2 <- pval2*length(pval2)/ranks2

for (r in length(pval1):1) {
  print(p_m_over_k1[ranks1>=r])
}

pval_adj1 <- pval_adj2 <- c()

for (i in 1:length(pval1)) {
  
  # find the rank
  rank1_temp <- ranks1[i]
  rank2_temp <- ranks2[i]
  
  # get all the p_m_over_k that are greater or equal to this rank and get the min value
  pval_adj1 <- c(pval_adj1, min(1,min(p_m_over_k1[ranks1>=rank1_temp])))
  pval_adj2 <- c(pval_adj2, min(1,min(p_m_over_k2[ranks2>=rank2_temp])))
}

### (values for) Table A.5 ###
cbind(label1, label2_1, pval1, ranks1, p_m_over_k1, pval_adj1) 

### (values for) Table A.6 ###
cbind(label1, label2_2, pval2, ranks2, p_m_over_k2, pval_adj2) 
# ####



# G. Attribute non-attendance (ANA) analysis (Table 6, 7, A.10, A.11, A.12) ---- 

eclc.mods <- list()

for(s in 1:8){
  
  #### b. ECLC-2 model (treatments) ______________________________________________
  #### 1. LOAD LIBRARY AND DEFINE CORE SETTINGS  ####
  ### Load Apollo library
  library(apollo)
  
  ### Initialise code
  apollo_initialise()
  
  ### Set core controls
  apollo_control = list(
    modelName  ="ECLC2-MXL",
    modelDescr ="ECLC model (MXL)",
    indivID    ="newid",
    nCores     = 10,
    mixing     = TRUE,
    panelData  = TRUE)
  # ####
  
  #### 2. LOAD DATA AND APPLY ANY TRANSFORMATIONS                     ####
  database <- read.csv("cost-vector-data.csv")
  database <- database[database$split==s, ]
  database$cost1 <- database$cost1/100
  database$cost2 <- database$cost2/100
  # ####
  
  #### 3. DEFINE MODEL PARAMETERS                                     ####

  if(s==1){
    # treatment 1
    apollo_beta <- c(-2.81194559506437, -26.2211189920332, 0.21418073174516, 0.610496620218221, 1.02308497704492, 0.555601009762126, 0.109281935274512, 1.47555090064097, 1.61738111169667, -0.0696502965206482, 0.374053139699254, 1.62906863076206, 1.68070558286744, -0.309801179680934, 2.4818175086198, 112.326738927756, 0.616226309524524, 0.00174325204151791, 0.506426636430138, 0.527135310237167, 1.15724945819635, 0.349324870831181, 0.46830339583528, 0.378098432492215, -0.819068304774814, -0.413479198462828, 0.0456718491402385, -1.70674841684687, 0.5)
  } else if(s==2){
    # treatment 2
    apollo_beta <- c(-2.61969333079792, -19.4672130151191, 0.371706707959431, 0.583984007885881, 0.83365184846754, 0.529567538869267, 0.0102396312112454, 1.39285334869965, 1.44346524564137, -0.151634228475704, 0.684591531405715, 1.78036411353003, 1.60910057353275, 0.0578946254682655, 2.39946974290469, 24.7735646454047, 0.420550898298913, -0.305828285556587, 0.140001860266222, -0.367398835555337, 0.774359584871565, 0.579124946942206, 0.578019676571201, -0.461780370620637, -0.148194680225686, 0.289625172466096, -0.487218171410803, -1.26674754133704, 0.5)
  } else if(s==3){  
    # treatment 3
    apollo_beta <- c(-2.52576941557431, -5.77300528126745, 0.57411854950652, 0.667039510644443, 1.00884536949696, 0.65391444638761, 0.389123939156325, 1.71087708165349, 1.4771016317228, -0.0679289609711944, 0.480650250324705, 1.73567446327696, 1.56739291487821, -0.788929249268449, 2.03202658134006, 16.6192450604111, 0.167463017466672, -0.462386721644203, -0.567909158755639, 0.53658062355908, -0.599229644240143, -0.459707506033695, 0.307672153218699, -0.394365062342123, -0.500904915417458, -0.892398124993979, 0.752329821510818, -1.75464116204381, 0.5)
  } else if(s==4){  
    # treatment 4
    apollo_beta <- c(-1.93305231713691, -10.6047172228997, -0.0175527000645109, 0.105162052423156, 0.554947530557239, 0.583379230794308, 0.52083670101741, 1.90011264773406, 1.1527923228688, -0.0931795155997277, 0.258695851659535, 1.37139841228807, 1.15694106971891, -1.07920705076379, 1.07647666378172, 11.9162402150509, -0.111387249886598, -0.382346008793829, -0.322042856828081, -0.234885851710937, 1.26146278241049, 0.156035267969607, -0.362623849334652, 0.211926199591096, -0.402605083460067, 0.123575179438341, -0.474500417767689, -2.13164093707177, 0.5)
  } else if(s==5){  
    # treatment 5
  apollo_beta <- c(-1.78837427118707, 0.114904123399546, 0.514353239446584, 0.89823649352916, 1.32070724289545, 0.581156229729382, -0.102859184013794, 1.34112447431278, 1.54948665385411, -0.253719311214405, 0.510251319347051, 1.28704724079716, 1.18947160571847, 1.5547032166921, 3.4174992103684, 1.41437190859366, -0.0465379207290801, -0.333242099182409, -0.782667576734176, -0.615100800064724, -1.38696579404395, 0.865428475519027, -0.222867232732934, -0.533903878191064, 0.311288876523567, -0.48015140255844, 0.837153663804297, -0.94113160312129, -2)
  } else if(s==6){  
    # treatment 6 
  apollo_beta <- c(-2.25344463421521, -0.0132956632093635, 0.446527283990117, 0.311420057718552, 0.711410140549512, 0.379393798843855, -0.362860303377584, 1.17222827052767, 1.27788354911574, -0.152584856844957, 0.0998157237521431, 1.09872877598465, 0.958149968883457, 1.22120392106549, 3.58859305614907, 1.4006069615754, 0.378243719652179, -0.93594854060565, 1.17064404036221, 0.436654410393568, 0.845452096735761, -0.489883254409479, -0.337200697902302, -0.400228748286837, 0.300141314286262, -0.51891448949601, 0.742455179539809, -1.1749733525448, -2) 
  } else if(s==7){  
    # treatment 7 
  apollo_beta <- c(-1.69613506738456, 0.55305211060566, 0.403531350837423, 0.694047079918518, 0.871702433640544, 0.456673668497237, -0.155544510968886, 0.9582423472983, 1.19124694564818, -0.0712627449344151, 0.335996467257638, 1.0248200764907, 0.928631321919078, 0.782502644129836, 1.52325078930647, 0.0308046589511, 0.0188200118002813, 0.554237050964608, -0.919407900412768, -0.0651315830869871, 0.454204523939422, 0.121167917780465, -0.34445105720345, -0.539607293776042, 0.0731556451687819, -0.236436094712658, 1.05235852024361, -1.88548066293832, -2.62663115814196)
  } else if(s==8){  
    # treatment 8 
  apollo_beta <- c(-2.09071691611866, 0.786695569298683, 0.66449327207786, 0.747981925502456, 0.60688915505103, 0.464091975045312, -0.134454361697733, 1.1385576201873, 1.0211865525014, -0.28094062687651, -0.162931402732858, 0.418183821053953, 0.753129175203865, 0.741976242065147, 1.36158411488601, -2.73352358934986, -0.911822823840043, 0.533475103496304, 1.07120405625837, 0.295278064731435, -0.248302409765031, 0.116933243001126, 0.0315916499719513, -0.178665497738286, 0.745868194505312, -0.572589493796532, -0.740651886966298, -2.00768185343074, -2.05586818342395) 
  }  
  
  names(apollo_beta) <- c("mu_none_a","mu_none_b",
                          "mu_clar1","mu_clar2","mu_clar3","mu_fish","mu_bio1","mu_bio2","mu_bio3","mu_coast","mu_lit1","mu_lit2","mu_lit3","mu_cost",
                          "sigma_none_a","sigma_none_b",
                          "sigma_clar1","sigma_clar2","sigma_clar3","sigma_fish","sigma_bio1","sigma_bio2","sigma_bio3","sigma_coast","sigma_lit1","sigma_lit2","sigma_lit3","sigma_cost",
                          "delta_b")
  apollo_fixed = c()   
  # ####
  
  #### 4.A. DEFINE RANDOM COMPONENTS                                    ####
  ### Set parameters for generating draws
  apollo_draws = list(
    interDrawsType = "sobol",
    interNDraws    = 1000,
    interNormDraws = c("draws_none_a", "draws_none_b", 
                       "draws_clar_1", "draws_clar_2", "draws_clar_3", "draws_fish",
                       "draws_bio_1", "draws_bio_2", "draws_bio_3", "draws_coast",
                       "draws_lit_1", "draws_lit_2", "draws_lit_3", "draws_cost")
  )
  
  ### Create random parameters
  apollo_randCoeff = function(apollo_beta, apollo_inputs){
    randcoeff = list()
    
    randcoeff[["b_cost_a"]]  = -exp(mu_cost + sigma_cost*draws_cost)
    randcoeff[["b_none_a"]]  = mu_none_a  + sigma_none_a*draws_none_a 
    randcoeff[["b_none_b"]]  = mu_none_b  + sigma_none_b*draws_none_b 
    randcoeff[["b_clar1_a"]] = mu_clar1 + sigma_clar1*draws_clar_1 
    randcoeff[["b_clar2_a"]] = mu_clar2 + sigma_clar2*draws_clar_2 
    randcoeff[["b_clar3_a"]] = mu_clar3 + sigma_clar3*draws_clar_3 
    randcoeff[["b_fish_a"]]  = mu_fish  + sigma_fish*draws_fish    
    randcoeff[["b_bio1_a"]]  = mu_bio1  + sigma_bio1*draws_bio_1   
    randcoeff[["b_bio2_a"]]  = mu_bio2  + sigma_bio2*draws_bio_2   
    randcoeff[["b_bio3_a"]]  = mu_bio3  + sigma_bio3*draws_bio_3   
    randcoeff[["b_coast_a"]] = mu_coast + sigma_coast*draws_coast  
    randcoeff[["b_lit1_a"]]  = mu_lit1  + sigma_lit1*draws_lit_1   
    randcoeff[["b_lit2_a"]]  = mu_lit2  + sigma_lit2*draws_lit_2   
    randcoeff[["b_lit3_a"]]  = mu_lit3  + sigma_lit3*draws_lit_3   
    
    return(randcoeff)
  }
  # ####
  
  #### 4.B. DEFINE LATENT CLASS COMPONENTS                              ####
  apollo_lcPars=function(apollo_beta, apollo_inputs){
    lcpars = list()
    
    lcpars[["b_none"]]  = list(b_none_a ,b_none_b )
    lcpars[["b_clar1"]] = list(b_clar1_a,b_clar1_a)
    lcpars[["b_clar2"]] = list(b_clar2_a,b_clar2_a)
    lcpars[["b_clar3"]] = list(b_clar3_a,b_clar3_a)
    lcpars[["b_fish"]]  = list(b_fish_a ,b_fish_a )
    lcpars[["b_bio1"]]  = list(b_bio1_a ,b_bio1_a )
    lcpars[["b_bio2"]]  = list(b_bio2_a ,b_bio2_a )
    lcpars[["b_bio3"]]  = list(b_bio3_a ,b_bio3_a )
    lcpars[["b_coast"]] = list(b_coast_a,b_coast_a)
    lcpars[["b_lit1"]]  = list(b_lit1_a ,b_lit1_a )
    lcpars[["b_lit2"]]  = list(b_lit2_a ,b_lit2_a )
    lcpars[["b_lit3"]]  = list(b_lit3_a ,b_lit3_a )
    lcpars[["b_cost"]]  = list(b_cost_a ,0        ) 
    
    V=list()
    V[["class_a"]] = 0
    V[["class_b"]] = delta_b # + delta_b_2*t2 + delta_b_3*t3 + delta_b_4*t4
    
    ### Settings for class allocation models
    classAlloc_settings = list(
      classes = c(class_a=1, class_b=2), 
      utilities = V
    )
    
    lcpars[["pi_values"]] = apollo_classAlloc(classAlloc_settings)
    
    return(lcpars)
  }
  # ####
  
  #### 5. GROUP AND VALIDATE INPUTS                                   ####
  apollo_inputs = apollo_validateInputs()
  # ####
  
  #### 6. DEFINE MODEL AND LIKELIHOOD FUNCTION                        ####
  apollo_probabilities=function(apollo_beta, apollo_inputs, functionality="estimate"){
    
    ### Attach inputs and detach after function exit
    apollo_attach(apollo_beta, apollo_inputs)
    on.exit(apollo_detach(apollo_beta, apollo_inputs))
    
    ### Create list of probabilities P
    P = list()
    
    ### Loop over classes
    for(s in 1:2){
      
      ### Define settings for MNL model component
      mnl_settings = list(
        alternatives = c(sq=1, alt2=2, alt3=3),
        avail        = 1,
        choiceVar    = choi
      )
      
      ### List of utilities: these must use the same names as in mnl_settings, order is irrelevant
      V = list()
      V[['sq']]   = b_cost[[s]]*cost0 + b_none[[s]] + b_clar1[[s]]*clar1_0 + b_clar2[[s]]*clar2_0 + b_clar3[[s]]*clar3_0 + b_fish[[s]]*fish0 + b_bio1[[s]]*bio1_0 + b_bio2[[s]]*bio2_0 + b_bio3[[s]]*bio3_0 + b_coast[[s]]*coast0 + b_lit1[[s]]*lit1_0 + b_lit2[[s]]*lit2_0 + b_lit3[[s]]*lit3_0
      V[['alt2']] = b_cost[[s]]*cost1               + b_clar1[[s]]*clar1_1 + b_clar2[[s]]*clar2_1 + b_clar3[[s]]*clar3_1 + b_fish[[s]]*fish1 + b_bio1[[s]]*bio1_1 + b_bio2[[s]]*bio2_1 + b_bio3[[s]]*bio3_1 + b_coast[[s]]*coast1 + b_lit1[[s]]*lit1_1 + b_lit2[[s]]*lit2_1 + b_lit3[[s]]*lit3_1
      V[['alt3']] = b_cost[[s]]*cost2               + b_clar1[[s]]*clar1_2 + b_clar2[[s]]*clar2_2 + b_clar3[[s]]*clar3_2 + b_fish[[s]]*fish2 + b_bio1[[s]]*bio1_2 + b_bio2[[s]]*bio2_2 + b_bio3[[s]]*bio3_2 + b_coast[[s]]*coast2 + b_lit1[[s]]*lit1_2 + b_lit2[[s]]*lit2_2 + b_lit3[[s]]*lit3_2
      
      mnl_settings$utilities = V
      mnl_settings$componentName = paste0("Class_",s)
      
      ### Compute within-class choice probabilities using MNL model
      P[[s]] = apollo_mnl(mnl_settings, functionality)
      
      ### Take product across observation for same individual
      P[[s]] = apollo_panelProd(P[[s]], apollo_inputs, functionality)
      
      ### Average across inter-individual draws within classes
      P[[s]] = apollo_avgInterDraws(P[[s]], apollo_inputs, functionality)
    }
    
    ### Compute latent class model probabilities
    lc_settings   = list(inClassProb = P, classProb=pi_values)
    P[["model"]] = apollo_lc(lc_settings, apollo_inputs, functionality)
    
    ### Prepare and return outputs of function
    P = apollo_prepareProb(P, apollo_inputs, functionality)
    return(P)
  }  
  # ####
  
  #### 7. MODEL ESTIMATION                                            ####
  estimate_settings = list(estimationRoutine="bgw", hessianRoutine="numDeriv")
  eclc.mods[[s]] = apollo_estimate(apollo_beta, apollo_fixed, apollo_probabilities, apollo_inputs, estimate_settings)
  
  # #####
}

### (values for) Table A.10 ###
apollo_modelOutput(eclc.mods[[1]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[2]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[3]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[4]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))

### (values for) Table A.11 ###
apollo_modelOutput(eclc.mods[[5]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[6]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[7]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[8]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))


# ALTERNATIVELY: Read in saved model estimates

eclc.mods <- readRDS("ECLC2-MXL-AllTreat.rds")

### (values for) Table A.10 ###
apollo_modelOutput(eclc.mods[[1]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[2]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[3]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[4]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))

### (values for) Table A.11 ###
apollo_modelOutput(eclc.mods[[5]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[6]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[7]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
apollo_modelOutput(eclc.mods[[8]], modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))





# Calculate marginal WTP and average posterior class membership probabilities for each model
nsim = 1000
treat.wtp.m <- treat.wtp.md <- list()
prob.class2 = matrix(0, nrow=8, ncol=3)

# Extract random draws from Apollo code above
emp_draws = emp_betas = list()
emp_draws[[1]]  = apollo_inputs$draws$draws_none_a[seq(1, nrow(apollo_inputs$draws$draws_none_a), 8),]
emp_draws[[2]]  = apollo_inputs$draws$draws_none_b[seq(1, nrow(apollo_inputs$draws$draws_none_b), 8),]
emp_draws[[3]]  = apollo_inputs$draws$draws_clar_1[seq(1, nrow(apollo_inputs$draws$draws_clar_1), 8),]
emp_draws[[4]]  = apollo_inputs$draws$draws_clar_2[seq(1, nrow(apollo_inputs$draws$draws_clar_2), 8),]
emp_draws[[5]]  = apollo_inputs$draws$draws_clar_3[seq(1, nrow(apollo_inputs$draws$draws_clar_3), 8),] 
emp_draws[[6]]  = apollo_inputs$draws$draws_fish[seq(1, nrow(apollo_inputs$draws$draws_fish), 8),]
emp_draws[[7]]  = apollo_inputs$draws$draws_bio_1[seq(1, nrow(apollo_inputs$draws$draws_bio_1), 8),]
emp_draws[[8]]  = apollo_inputs$draws$draws_bio_2[seq(1, nrow(apollo_inputs$draws$draws_bio_2), 8),]
emp_draws[[9]]  = apollo_inputs$draws$draws_bio_3[seq(1, nrow(apollo_inputs$draws$draws_bio_3), 8),]
emp_draws[[10]] = apollo_inputs$draws$draws_coast[seq(1, nrow(apollo_inputs$draws$draws_coast), 8),]
emp_draws[[11]] = apollo_inputs$draws$draws_lit_1[seq(1, nrow(apollo_inputs$draws$draws_lit_1), 8),]
emp_draws[[12]] = apollo_inputs$draws$draws_lit_2[seq(1, nrow(apollo_inputs$draws$draws_lit_2), 8),]
emp_draws[[13]] = apollo_inputs$draws$draws_lit_3[seq(1, nrow(apollo_inputs$draws$draws_lit_3), 8),]
emp_draws[[14]] = apollo_inputs$draws$draws_cost[seq(1, nrow(apollo_inputs$draws$draws_cost), 8),]

for(t in 1:8){
  model <- eclc.mods[[t]]
  apollo_modelOutput(model, modelOutput_settings = list(printPVal=TRUE, printClassical=FALSE))
  betas = model$estimate[-c(1:2,15:16)]
  prob_cl1 = 1 / (1 + exp(model$estimate["delta_b"]))
  draws_temp = mvrnorm(nsim, mu=model$estimate, Sigma=model$robvarcov)
  treat.wtp.m_temp <- treat.wtp.md_temp <- matrix(0, nrow=11, ncol=3)
  
  # obtain unconditional parameter distributions
  for(i in 1:11){
    emp_betas[[i]] = betas[i] + betas[i+12] * emp_draws[[i+2]]
  }
  emp_betas[[12]] = -exp(betas[12] + betas[24] * emp_draws[[14]])    # cost
  
  # compute WTP
  for(i in 1:11){
    treat.wtp.m_temp[i,1]  = mean(-100*emp_betas[[i]] / emp_betas[[12]])
    treat.wtp.md_temp[i,1] = median(-100*emp_betas[[i]] / emp_betas[[12]])
  }
  
  # Simulate WTP confidence intervals
  emp_betas_temp <- list()
  wtp.m_temp <- wtp.md_temp <- matrix(0, nrow=nsim, ncol=11)
  for(s in 1:nsim){
    betas_temp = draws_temp[s,-c(1:2,15:16)]
    for(l in 1:11){
      emp_betas_temp[[l]] = betas_temp[l] + betas_temp[l+12] * emp_draws[[l+2]]
    }
    emp_betas_temp[[12]] = -exp(betas_temp[12] + betas_temp[24] * emp_draws[[14]])    # cost
    
    # compute wtp
    for(l in 1:11){
      wtp.m_temp[s,l]  = mean(-100*emp_betas_temp[[l]] / emp_betas_temp[[12]])
      wtp.md_temp[s,l] = median(-100*emp_betas_temp[[l]] / emp_betas_temp[[12]])
    }
  }
  for(l in 1:11){
    treat.wtp.m_temp[l,2:3] = quantile(wtp.m_temp[,l], probs=c(0.025,0.975))
    treat.wtp.md_temp[l,2:3] = quantile(wtp.md_temp[,l], probs=c(0.025,0.975))
  }
  # store treatment-specific WTP matrices
  treat.wtp.m[[t]] = treat.wtp.m_temp
  treat.wtp.md[[t]] = treat.wtp.md_temp
  
  rownames(treat.wtp.m[[t]]) = rownames(treat.wtp.md[[t]]) = c("clar1","clar2","clar3","fish","bio1","bio2","bio3","coast","lit1","lit2","lit3")
  colnames(treat.wtp.m[[t]])  = c("Mean","lb","ub")
  colnames(treat.wtp.md[[t]]) = c("Median","lb","ub")
  
  # Average class membership probability
  prob.class2[t,1] = 1 - prob_cl1
  prob_cl2_dist = exp(draws_temp[,29]) / (1 + exp(draws_temp[,29]))
  prob.class2[t,2:3] = quantile(prob_cl2_dist, probs=c(0.025,0.975))
}

colnames(prob.class2) = c("mean.cl2","lb","ub")
rownames(prob.class2) = c("T1","T2","T3","T4","T5","T6","T7","T8")

eclc.median.wtp <- matrix(0, nrow=11, ncol=8)
for(t in 1:8){
  eclc.median.wtp[,t] = treat.wtp.md[[t]][,1]
}
colnames(eclc.median.wtp) <- c("T1","T2","T3","T4","T5","T6","T7","T8")
rownames(eclc.median.wtp) <- c("clar1","clar2","clar3","fish","bio1","bio2","bio3","coast","lit1","lit2","lit3")


### Table A.12 ###
round(eclc.median.wtp, 0)

### Table 6 ###
round(prob.class2, 3)

# ####



# H. Stated attribute non-attendance (ANA) (Table 7) ---------------------------

database <- read.csv("cost-vector-data.csv")
data <- database[database$task==1,]

# Percentages stating "Never" or "Always" per treatment
stated.ana = matrix(0, ncol=8, nrow=2)

for(s in 1:8){
  stated.ana[1,s] = (table(data$v_390[data$split==s]==1) / length(data$v_390[data$split==s]))[2]  # never attended cost (1)
  stated.ana[2,s] = (table(data$v_390[data$split==s]==5) / length(data$v_390[data$split==s]))[2]  # always attended cost (5) 
}
colnames(stated.ana) <- c("T1","T2","T3","T4","T5","T6","T7","T8")
rownames(stated.ana) <- c("never","always")
round(stated.ana*100, 2)

chi2.never = chi2.always = matrix(0, nrow=4, ncol=1)

for(s in 1:4){
  to.test = rbind(table(data$v_390[data$split==s]==1), table(data$v_390[data$split==(4+s)]==1))
  test = chisq.test(to.test)
  chi2.never[s,] = test$p.value
  
  to.test = rbind(table(data$v_390[data$split==s]==5), table(data$v_390[data$split==(4+s)]==5))
  test = chisq.test(to.test)
  chi2.always[s,] = test$p.value
}
rownames(chi2.never) = rownames(chi2.always) = c("T1 vs T5","T2 vs T6","T3 vs T7","T4 vs T8")
colnames(chi2.never) = colnames(chi2.always) = "p.value"

### values for Table 7 ###
round(stated.ana*100, 2)
round(chi2.never, 3)
round(chi2.always, 3)

# ####



# I. Realism and consequentiality items (Table 8, A.7) ------------------------- 

database <- read.csv("cost-vector-data.csv")
data <- database[database$task==1,]

item.a <- item.b <- item.c <- item.d <- item.e <- item.f <- list()

for(i in 1:8){
  item.a[[i]] = table(data$v_627[data$split==i])
  item.b[[i]] = table(data$v_628[data$split==i])
  item.c[[i]] = table(data$v_631[data$split==i])
  item.d[[i]] = table(data$v_622[data$split==i])
  item.e[[i]] = table(data$v_624[data$split==i])
  item.f[[i]] = table(data$v_629[data$split==i])
}


### (values for) Table 8 ###

fishers_exact_pvalues <- list()
set.seed(12354)

# A. I believe that the money from the Baltic Sea tax would actually be spent on the protection of the Baltic Sea.
to.test = rbind(item.a[[1]],item.a[[5]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[1]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.a[[2]],item.a[[6]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[2]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.a[[3]],item.a[[7]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[3]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.a[[4]],item.a[[8]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[4]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

# B. I believe that the programmes announced will actually contribute to a better state of the Baltic Sea.
to.test = rbind(item.b[[1]],item.b[[5]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[5]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.b[[2]],item.b[[6]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[6]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.b[[3]],item.b[[7]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[7]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.b[[4]],item.b[[8]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[8]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value


# C. I have been thinking about what I could no longer afford for myself and/or my and/or my family if I spend money on for the protection of the Baltic Sea.
to.test = rbind(item.c[[1]],item.c[[5]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[9]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.c[[2]],item.c[[6]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[10]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.c[[3]],item.c[[7]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[11]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.c[[4]],item.c[[8]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[12]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value


# D. I would be willing to pay any price so a good environmental status in the Baltic Sea is achieved; and I can afford the payments.
to.test = rbind(item.d[[1]],item.d[[5]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[13]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.d[[2]],item.d[[6]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[14]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.d[[3]],item.d[[7]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[15]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.d[[4]],item.d[[8]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[16]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value


# E. My decisions will influence whether and to what extent the programmes are implemented.
to.test = rbind(item.e[[1]],item.e[[5]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[17]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.e[[2]],item.e[[6]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[18]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.e[[3]],item.e[[7]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[19]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.e[[4]],item.e[[8]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[20]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value


# F. The survey is to find out how the public's attitude towards the programmes to protect the Baltic Sea; this does not imply a commitment to pay for them.
to.test = rbind(item.f[[1]],item.f[[5]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[21]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.f[[2]],item.f[[6]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[22]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.f[[3]],item.f[[7]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[23]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value

to.test = rbind(item.f[[4]],item.f[[8]])
t(to.test)
round(t(rbind(to.test[1,] / sum(to.test[1,]) , to.test[2,] / sum(to.test[2,]))), 2)
fisher.test(to.test, simulate.p.value=TRUE)$p.value
fishers_exact_pvalues[[24]] = fisher.test(to.test, simulate.p.value=TRUE)$p.value


# p-values of Fisher's exact tests (above)
pvalues <- unlist(fishers_exact_pvalues)

label1 <- rep(c("statement A","statement B","statement C","statement D","statement E","statement F"), each=4) 
label2 <- rep(c("T1 vs. T5","T2 vs. T6","T3 vs. T7","T4 vs. T8"), 6)

ranks<-rank(pvalues, ties.method = "last")
p_m_over_k<-pvalues*length(pvalues)/ranks
cbind(pvalues, ranks, p_m_over_k) 

for (r in length(pvalues):1) {
  print(p_m_over_k[ranks>=r])
}

pvalues_adj<-c()

for (i in 1:length(pvalues)) {
  
  # find the rank
  tmp_rank<-ranks[i]
  
  # get all the p_m_over_k that are greater or equal to this rank and get the min value
  pvalues_adj<-c(pvalues_adj, min(1,min(p_m_over_k[ranks>=tmp_rank])))
}
print(pvalues_adj)

### Table A.7 ###
round(cbind(pvalues, ranks, p_m_over_k, pvalues_adj), 4)

# ####

